home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / test2.i < prev    next >
Text File  |  1995-07-26  |  12KB  |  369 lines
  1. /*
  2.    TEST2.I
  3.    Highly vectorizing physics problem for timing Yorick.
  4.  
  5.    $Id: test2.i,v 1.1 1993/08/27 18:50:06 munro Exp $
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. if(!is_void(is_a_mac)) {
  11.   /* reduce the problem size on the Mac so I don't have to increase 
  12.      the default memory. */
  13.   big_test= 0;
  14. } else {
  15.   big_test= 1;
  16. }
  17.  
  18. func test2(npass)
  19. /* DOCUMENT test2
  20.          or test2, npass
  21.      Given a slab divided into zones by parallel planes, and given a
  22.      set of photon group boundary energies, compute the contribution of
  23.      each zone to the radiation flux emerging from one surface of the
  24.      slab.  If NPASS is given, the calculation is repeated that many
  25.      times.  The zoning, photon group structure, opacities, and source
  26.      functions are all computed arbitrarily, but the number of zones
  27.      and groups are taken to be representative of a typical 1-D
  28.      radiation transport calculation.  */
  29. {
  30.   if (is_void(npass) || npass<=0) npass= 1;
  31.  
  32.   now= split= exponential= array(0.0, 3);
  33.   timer, now;
  34.   n= npass;
  35.   while (n--) esc= escout(zt, akap, srcfun);
  36.   timer, now, split;
  37.   timer_print, "Time per pass", split/npass, "Total time", split,\
  38.     "Computing exponentials", exponential;
  39.  
  40.   if(big_test > 0) {
  41.     if (esc(0,mxx)!=29 || !approx_eq(esc(0,max),0.001875213417) ||
  42.         !approx_eq(esc(0,1),0.001694423767) ||
  43.         !approx_eq(esc(0,0),5.434635527e-05) ||
  44.         esc(1,mxx)!=39 || !approx_eq(esc(1,max),0.000440460566))
  45.       write, "***WARNING*** values returned by escout are not correct";
  46.  
  47.     if (flxout(0,mxx)!=34 || !approx_eq(flxout(0,max),0.06615472064) ||
  48.         !approx_eq(flxout(0,1),0.003187516911) ||
  49.         !approx_eq(flxout(0,0),0.005280842058) ||
  50.         flxout(1,mxx)!=29 || !approx_eq(flxout(1,max),0.001805157164))
  51.       write, "***WARNING*** values of flxout are not correct";
  52.  
  53.     if (!approx_eq(min(tau(:-1,)),6.982090961e-05) ||
  54.         !approx_eq(max(tau),45.80160946))
  55.       write, "***WARNING*** values of tau are not correct";
  56.   }
  57. }
  58.  
  59. func approx_eq(x, y)
  60. {
  61.   return (abs(x-y)/(abs(x+y)+1.e-35))<1.e-6;
  62. }
  63.  
  64. #if 0
  65. CPU seconds per pass:
  66.             IDL    Yorick     Basis       DAP     FORTRAN(-O)
  67. HP730        -      0.60       2.04      1.84       0.27  (0.60 -g)
  68. Solbourne  2.81     1.90       6.02      5.90       1.00
  69.     (varies by ~10%)
  70.  
  71. Using forum (Solbourne) on 8/Dec/92:
  72.      forum[18] time /usr/local/lib/idl/bin.sunos.4.1.sun4/idl
  73.      IDL. Version 2.3.0 (sunos sparc).
  74.      Copyright 1989-1992, Research Systems, Inc.
  75.      All rights reserved.  Unauthorized reproduction prohibited.
  76.      Site: 1491.
  77.      Licensed for use by: LLNL - X Division
  78.      IDL> .rnew /home/miggle/munro/Yorick/include/test2
  79.      % Compiled module: ESCOUT.
  80.      % Compiled module: APPROX_EQ.
  81.      % Compiled module: TEST2.
  82.      % Compiled module: BNU.
  83.      % Compiled module: OPACSET.
  84.      % Compiled module: SPAN.
  85.      % Compiled module: SPANL.
  86.      % Compiled module: $MAIN$.
  87.      IDL> test2,1
  88.      IDL> exit
  89.      3.2u 0.4s 0:12 29% 0+2496k 1+2io 1pf+0w
  90.      forum[19] time time /usr/local/lib/idl/bin.sunos.4.1.sun4/idl
  91.      IDL. Version 2.3.0 (sunos sparc).
  92.      Copyright 1989-1992, Research Systems, Inc.
  93.      All rights reserved.  Unauthorized reproduction prohibited.
  94.      Site: 1491.
  95.      Licensed for use by: LLNL - X Division
  96.      IDL> .rnew /home/miggle/munro/Yorick/include/test2
  97.      % Compiled module: ESCOUT.
  98.      % Compiled module: APPROX_EQ.
  99.      % Compiled module: TEST2.
  100.      % Compiled module: BNU.
  101.      % Compiled module: OPACSET.
  102.      % Compiled module: SPAN.
  103.      % Compiled module: SPANL.
  104.      % Compiled module: $MAIN$.
  105.      IDL> test2,11
  106.      IDL> exit
  107.      31.3u 0.9s 0:58 54% 0+3544k 0+0io 0pf+0w
  108.  
  109.  
  110.      forum[20] time yorick
  111.       Yorick ready.  For help type 'help'
  112.      > #include "/home/miggle/munro/Yorick/include/test2.i"
  113.      > test2,1
  114.             Timing Category     CPU sec  System sec    Wall sec
  115.               Time per pass       1.840       0.300       2.270
  116.              Total time       1.840       0.300       2.270
  117.          Computing exponentials       0.600       0.120       0.760
  118.       -----Total Elapsed Times-----       2.510       0.490      27.500
  119.      > test2,10
  120.             Timing Category     CPU sec  System sec    Wall sec
  121.               Time per pass       1.932       0.050       1.986
  122.              Total time      19.320       0.500      19.860
  123.          Computing exponentials       6.360       0.080       6.440
  124.       -----Total Elapsed Times-----      21.860       1.020      53.080
  125.      > quit
  126.      21.9u 1.0s 0:59 38% 0+3560k 31+0io 85pf+0w
  127.  
  128.  
  129.      forum[21] basis
  130.      Basis    (basis, Version 921125)
  131.      Run at 13:02:49 on 12/08/92 on the forum    machine, suffix 10797x
  132.      Initializing Basis System
  133.      Basis 7.0
  134.      Initializing EZCURVE/NCAR Graphics
  135.      Ezcurve/NCAR 2.0
  136.      Basis> echo=0
  137.      Basis> read /home/miggle/munro/Yorick/include/test2.bas
  138.      End of input file /home/miggle/munro/Yorick/include/test2.bas
  139.      Resuming input from TERMINAL
  140.      Basis> test2(1)
  141.  
  142.  
  143.     CPU (sec)   SYS (sec)
  144.         6.017       0.300
  145.      Basis> test2(10)
  146.  
  147.     CPU (sec)   SYS (sec)
  148.        60.233       0.017
  149.      Basis> end
  150.  
  151.     CPU (sec)   SYS (sec)
  152.        68.617       0.917
  153.      68.6u 0.9s 2:07 54% 0+6144k 19+1io 213pf+0w
  154.  
  155.  
  156.      forum[22] dap
  157.      DAP> read /home/miggle/munro/Yorick/include/test2.bas
  158.      DAP> test2(1)
  159.  
  160.     CPU (sec)   SYS (sec)
  161.         5.883       0.400
  162.      DAP> test2(10)
  163.  
  164.     CPU (sec)   SYS (sec)
  165.        59.533       0.000
  166.      DAP> end
  167.      69.1u 1.5s 1:47 65% 0+6488k 20+1io 218pf+0w
  168.  
  169.  
  170. Yorick results (test2.bas) on tonto (HP730) 21:53 4/Dec/92
  171. top showed SIZE/RES  764K/244K  at prompt
  172.                     3636K/3052K after test2,1
  173.                     4628K/4012K after test2,20
  174.      tonto[9] yorick
  175.       Yorick ready.  For help type 'help'
  176.      > #include "/home/miggle/munro/Yorick/include/test2.i"
  177.      > test2,1
  178.             Timing Category     CPU sec  System sec    Wall sec
  179.               Time per pass       0.550       0.100       0.670
  180.              Total time       0.550       0.100       0.670
  181.          Computing exponentials       0.240       0.030       0.270
  182.       -----Total Elapsed Times-----       0.930       0.220      34.950
  183.      > test2,20
  184.             Timing Category     CPU sec  System sec    Wall sec
  185.               Time per pass       0.599       0.001       0.906
  186.              Total time      11.980       0.020      18.110
  187.          Computing exponentials       5.030       0.000       7.270
  188.       -----Total Elapsed Times-----      12.920       0.240      83.120
  189.  
  190. Basis results (test2.bas) on tonto (HP730) 21:53 4/Dec/92
  191. top showed SIZE/RES 5356K/540K  at prompt
  192.                     8960K/4108K after test2(1)
  193.                     8960K/4112K after test2(20)
  194.      tonto[8] basis
  195.      Basis    (basis, Version 921125)
  196.      Run at 21:45:42 on 12/04/92 on the tonto    machine, suffix 27713x
  197.      Initializing Basis System
  198.      Basis 7.0
  199.      Initializing EZCURVE/NCAR Graphics
  200.      Ezcurve/NCAR 2.0
  201.      Basis> echo=0
  202.      Basis> read /home/miggle/munro/Yorick/include/test2.bas
  203.      End of input file /home/miggle/munro/Yorick/include/test2.bas
  204.      Resuming input from TERMINAL
  205.      Basis> test2(1)
  206.  
  207.     CPU (sec)   SYS (sec)
  208.         1.990        .120
  209.      Basis> test2(20)
  210.  
  211.     CPU (sec)   SYS (sec)
  212.        40.800        .000
  213.  
  214. DAP results (test2.bas) on tonto (HP730) 21:53 4/Dec/92
  215. top showed SIZE/RES 6796K/844K   at prompt
  216.                     10360K/4352K after test2(1)
  217.                     10360K/4352K after test2(20)
  218.      tonto[11] dap
  219.      DAP> echo=0
  220.      DAP> read /home/miggle/munro/Yorick/include/test2.bas
  221.      DAP> test2(1)
  222.  
  223.     CPU (sec)   SYS (sec)
  224.         1.840        .080
  225.      DAP> test2(20)
  226.  
  227.     CPU (sec)   SYS (sec)
  228.        36.830        .000
  229. #endif
  230.  
  231. /* This routine computes the optical depth through each zone at every
  232.    frequency and then uses that to compute the radiation emitted in
  233.    each zone that escapes from the problem, assuming plane parallel
  234.    geometry.
  235.    The returned result is an nzones-by-ngroups array of
  236.    power per unit photon energy per unit area.  */
  237. func escout(zt,      /* npoints zone boundary positions (cm) */
  238.         akap,    /* nzones-by-ngroups opacities (1/cm) */
  239.         srcfun)  /* nzones-by-ngroups source (specific intensity units) */
  240. {
  241.   extern flxout;     /* (output) nzones-by-ngroups outgoing fluxes */
  242.   extern dtau;       /* (output) nzones-by-ngroups optical depths (ODs) */
  243.   extern tau;        /* (output) npoints-by-ngroups cumulative ODs */
  244.  
  245.   extern mu, wmu;    /* Gauss-Legendre cos(theta) and weight arrays
  246.                 for integration over escape angles */
  247.  
  248.   /* compute tau, the optical depth to each zone along the zt-direction */
  249.   dtau= akap*zt(dif);
  250.   tau= array(0.0, dimsof(zt), dimsof(akap)(3));
  251.   tau(2:,)= dtau(psum,);
  252.   /* consider the outward going radiation, and thus use tau
  253.      measured from the right boundary */
  254.   tau= tau(0,)(-,) - tau;
  255.  
  256.   /* compute exf, the fraction of the srcfun which escapes from
  257.      each zone in each bin at each angle mu relative to the zt direction */
  258.   enow= array(0.0, 3);
  259.   timer, enow;  /* most of the actual work is computing these exponentials */
  260.   exf= exp(-tau(,,-)/mu(-,-,));
  261.   timer, enow, exponential;
  262.  
  263.   /* compute the escaping flux per unit surface area in each bin
  264.      contributed by each zone -- units are 10^17 W/kev/cm^2  */
  265.   /* Note:
  266.        dI/dtau = -I + S   (S is srcfun, which is Bnu) implies
  267.        I = integral of ( S * exp(tau) ) dtau from tau -infinity to tau 0
  268.        Hence, the contribution of a single zone to this integral is
  269.        S12 * (exp(tau1) - exp(tau2)), which explains the exf(dif,,) below.
  270.        The wmu are Gauss-Legendre integration weights, and the mu is the
  271.        cos(theta) to project the specific intensity onto the direction
  272.        normal to the surface (zt).  */
  273.   esfun= 2*pi*exf(dif,,)*srcfun(,,-)*(mu*wmu)(-,-,);
  274.  
  275.   /* Also compute what the total flux WOULD have been if each successive
  276.      zone were at the surface.  This is the same as the "one-sided"
  277.      flux directed outward across each zone boundary.  */
  278.   fuzz= 1.0e-10;
  279.   flxout= (esfun(psum,,)/(exf(2:,,)+fuzz))(,,sum);
  280.  
  281.   return esfun(,,sum);
  282. }
  283.  
  284. /* This function is used to set a dummy opacity to be used by the
  285.    transport calculation  */
  286. func opacset(tmp, rho, gav, gb)
  287. {
  288.   extern srcfun;  /* (output) nzones-by-ngroups Bnu, LTE source function */
  289.   extern akap;    /* (output) nzones-by-ngroups opacities (1/cm) */
  290.  
  291.   /* the opacity is proportional to the density to the rhopow
  292.      power and the temperature to the tempow power  */
  293.   rhopow= 2;
  294.   tempow= 3;
  295.   factr= (tmp/1.0)^tempow*(rho/1.0)^rhopow;
  296.  
  297.   /* set the constants in front of the terms for the two edges */
  298.   con0= 1.0e+4;
  299.   cona= 2.5e+2;
  300.   conb= 5.0e+0;
  301.   /* set the energies of the two edges */
  302.   edg0= 0.1;
  303.   edga= 0.5;
  304.   edgb= 2.0;
  305.  
  306.   /* set up arrays that are zero below the edges and one above them  */
  307.   vala= double(gav>edga);
  308.   valb= double(gav>edgb);
  309.  
  310.   /* frequency dependence is the same for all zones, so calculate it once */
  311.   frval= con0*(edg0/gav)^3+cona*vala*(edga/gav)^3+conb*valb*(edgb/gav)^3;
  312.  
  313.   /* set the opacity */
  314.   akap= factr(,-)*frval(-,);
  315.   /* the source function is always a blackbody */
  316.   srcfun= bnu(tmp, gb);
  317. }
  318.  
  319. /* Make a blackbody using the analytic formula --
  320.    the units are 10^17 W/kev/g/ster */
  321. func bnu(tmp, freqb)
  322. {
  323.   /* Note that tmp and freqb are one dimensional and the output
  324.      array should be 2D  */
  325.   /* compute the derivative using the exact black body,
  326.      not the average over the bin  */
  327.   exf= exp(-(freqb(zcen)(-,)/tmp));
  328.   return 0.0504*(freqb(zcen)^3)(-,)*exf/(1.0-exf);
  329. }
  330.  
  331. /* set up default params. etc. */
  332.  
  333. /* set the number of spatial zones and photon groups */
  334. if(big_test > 0) {
  335.   npoints= 100;
  336. } else {
  337.   npoints= 40;
  338. }
  339. if(big_test >= 0) {
  340.   ngroups= 64;
  341. } else {
  342.   ngroups= 32;
  343. }
  344. nzones= npoints-1;
  345.  
  346. /* set the zone boundary positions */
  347. zt= span(0.0, 0.0050, npoints);
  348.  
  349. /* set the frequency bins */
  350. gb= spanl(0.1, 4.0, ngroups+1);
  351. gav= gb(zcen);
  352.  
  353. /* set the density and temperature */
  354. rho= spanl(1.0, 1.0, nzones);
  355. tmp= spanl(1.0, 1.0, nzones);
  356.  
  357. /* set the opacity akap and source function srcfun */
  358. opacset, tmp, rho, gav, gb;
  359.  
  360. /* set up for a Gauss-Legendre angular quadrature */
  361. xmu= [0.1488743389, 0.4333953941, 0.6794095682, 0.8650633666, 0.9739065285];
  362. mu= -xmu(::-1);
  363. grow, mu, xmu;
  364. xmu= [0.2955242247, 0.2692667193, 0.2190863625, 0.1494513491, 0.0666713443];
  365. wmu= xmu(::-1);
  366. grow, wmu, xmu;
  367. /* correct the Gauss-Legendre abcissas-- interval is only 0 to 1 */
  368. mu= 0.5 + 0.5*mu;
  369.